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Abstract. - We study molecular diffusion in linear nanopores with different types of roughness 
in the so-called Knudsen regime. Knudsen diffusion represents the limiting case of molecular 
diffusion in pores, where mutual encounters of the molecules within the free pore space may be 
neglected and the time of flight between subsequent collisions with the pore walls significantly 
exceeds the interaction time between the pore wall and the molecules. We present an extension 
of a commonly used procedure to calculate transport diffusion coefficients. Our results show that 
using this extension, the coefficients of self- and transport diffusion in the Knudsen regime are 
equal for all regarded systems, which improves previous literature data. 
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Introduction. — Diffusion phenomena of gases in dis- 
ordered and porous media have been subject to intense re- 
search for several decades [1-4] with applications in hetero- 
geneous catalysis [5] , adsorption [6] and separation [7] . Re- 
cent progress in synthesizing nanostructured porous ma- 
terials [6,8] has provided essentially unlimited options for 
the generation of purpose-taylored pore architectures and 
there is an increasing demand for clarification of the main 
features of molecular transport in such systems [9, 10]. In 
this work, we concentrate on transport pores, that play 
an important role in bimodal porous materials [11] where 
they ensure fast molecular exchange between the micro- 
porous regions, in which the actual conversion and sepa- 
ration phenomena take place. In the transport pores, the 
so-called Knudsen diffusion dominates, where the inter- 
molecular collisions can be neglected and the molecules 
perform a series of free flights and change direction statis- 
tically after collisions with the pore walls. 

Experimentally, two kinds of diffusion problems can 
be considered, the so-called transport diffusion, where 
the particles diffuse in a non-equilibrium situation from 
one side of the system to the opposite side (here under 
the influence of a concentration gradient) and the self- 
(or tracer-) diffusion under equilibrium conditions. Both 
problems are described by the transport diffusion coeffi- 
cient D t and the self- (or tracer-) diffusion coefficient D s , 
respectively. D t is defined by Fick's 1st law as the propor- 
tionality constant between the current density j and the 
concentration gradient dc/dx, 



< x 2 (t) > of a random walker 
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while D s is defined by the moan square displacement 



after time t, where d is the dimension of the pore. Under 
the conditions of Knudsen diffusion, both diffusion coeffi- 
cients are expected to coincide [2]. Using the method of 
Evans et al [12], this has recently been verified numerically 
for smooth pores and for pores of different surface rough- 
ness [13] showing, however, slight numerical differences 
between D s and D t . Since, on the other side, there are 
also studies which do in fact consider the possibility of a 
difference between transport and self-diffusion coefficients 
in the Knudsen regime [15], even these differences have to 
be taken seriously. In [13] their occurence is attributed to 
the possibility of a non-constant density gradient of the 
gas particles inside the pores that disturbs the method 
of [12] slightly. Here, we investigate this problem in de- 
tail both for diffusion in the pores and for quite general 
random walks. We find that the deviations from the ex- 
pected gradient are strongest when the particles perform 
discrete jumps with many different velocities. The devia- 
tions disappear completely, if only one fixed jump length 
per time step is allowed and they are quite small (but mea- 
surable) for the Knudsen pores of Ref. [13] (see Fig. 1). 
We present an improvement of Evan's method (that we 
call "Enhanced f t Method" (EFM)) that accounts for the 
non-constant gradients and show that with EFM, D s = D t 
for all considered systems. 

Calculation of the transport diffusion. In the 

simulations of D t as well as in a typical diffusion exper- 
iment, a concentration gradient dc(x)/dx is applied with 
the concentrations c = cq = 1 on the left hand side and 
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Fig. 1: Realizations of the pore geometry generated by a gen- 
eralized random 3d-Koch curve of generation v. (a) Smooth 
pore (y = 0) and (b-c) rough pores of v = 1 and 2 of length 
L and width h (for details see [13]). For the simulations, the 
system is covered with a grid of lattice constant a — h/64. 

c = on the right hand side (x > L) of the pore, respec- 
tively. All particles start at x = and perform a random 
walk in d = 1 or a random trajectory between the system 
walls in d — 3, using Lambert's law of reflection [13]. They 
are absorbed when they hit the left or the right boundary 
(Dirichlct boundary conditions). After some relaxation 
time, this leads to a constant current density j, described 
by Eq. (1). 

Relaxation of a particle flow into a stationary state is 
very time-consuming, as the particle flux has to be moni- 
tored throughout the system. In the stationary state the 
particle concentration does not depend on the relaxation 
time t. In our simulations we consider a state as sta- 
tionary if the particle concentrations between two given 
time steps differ at maximum by a predetermined thresh- 
old value [14]. It is common practice to use the faster 
Evans' method [12] to derive j from the (transmission) 
probability f t that a particle starting at the left boundary 
will leave the pore through the right boundary, 

j = c ft <u x > (3) 

where < u x > is the mean velocity in x-direction. Com- 
bining Eqs. (1) and (3) yields 

A = -co/* < u x > (|) . (4) 

Usually, the concentration gradient is assumed to be con- 
stant and equal to 



Combining Eqs. (4) and (5) yields 

D t = ft<u x > L. (6) 

Accordingly, for calculating f t , N random trajectories are 
considered that start at x = and end when cither x = 
or x = L is reached. As we show in this paper, the prob- 
lem with Evans' method is that it only works if the con- 
centration gradient is well described by Eq. (5). It leads 
to spurious results, if deviations from a constant concen- 
tration gradient occur and in this case, Eq. (4) has to be 
taken as starting point. 
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Fig. 2: Concentration densities c(x/h) as calculated through 
relaxation of 10 6 particles into a stationary state (open sym- 
bols) and by the Enhanced f t Method (EFM, filled symbols), 
plotted versus x/hior smooth 3d pores of length L = 50/i. The 
results were generated from 10 6 runs. 

From probability density to concentration. To 

obtain the correct value of D t via Eq. (4), we need the con- 
centration c(x) within the pore and the associated concen- 
tration gradient dc(x)/dx. Figure 2 shows the results for 
c(x) as calculated by relaxing the density of gas particles 
inside smooth 3d pores of length L — 50h and width h 
into a stationary state (open symbols). 

For the simulation h was subdivided into smaller seg- 
ments of size a with h = 64a. Therefore, a is the lattice 
constant which for the roughest pores considered is the 
size of the smallest boundary structures. When the par- 
ticle hits a wall, the new direction is chosen randomly, 
according to Lambert's law of reflection. The path of the 
particle and hence the next collision is calculated by using 
basic mathematics within successive volume elements of 
size a 3 . The traveled distances are evaluated after all inte- 
ger numbers of time steps r, whereas for the path lengths, 
non-integer multiples of a are allowed and determined by 
linear interpolation. 

We can see that the concentration profile in Fig. 2 dif- 
fers only inside a small boundary region from a linear 
profile. Close to the left boundary, c(x) shows a small 
bump, whose relative position x m - ln /L approaches zero for 
increasing system size L. Therefore, we can still apply 
Evans' method for x ^> x m i n , where dc(x)/dx = const but 
unequal to —1/L. The deviation of the slope from — 1/L 
changes D t by several percent. Unfortunately, the rela- 
tive deviation between the correct dc(x)/dx and —1/L is 
unaffected by the system size, as the absolute deviation de- 
creases with the system size in the same way as dc(x)/dx. 
Hence, the problem cannot be solved by simply increasing 
the system size, but only by calculating the correct c(x). 

Therefore, if we want to apply Evans method, we must 
find a fast procedure to calculate the correct c{x) for large 
systems. Figure 3 shows the schematics of such a pro- 
cedure in Id and 2d, where c(x) is calculated from the 
trajectories used to calculate ft- The particle positions at 
every time step are indicated by dots. For better statis- 
tics we used 64 time steps in every time interval r that 
a particle with u = 1 needs to travel the distance of h. 
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Fig. 3: Schematics for generating the particle histogram and 
hence c(x) for (a) Id jumps with constant velocity and (b) 
particle trajectories in 2D. Every time step, the histograms 
(below) are incremented at the appropriate location. 



This has to be taken into account when normalizing the 
concentration c(x). The histogram over all particle po- 
sitions of the whole simulation is shown at the bottom 
of the figure. Hence, the histogram describes the time- 
averaged concentration c(x) . Using the ergodic hypothesis 
we identify c(x) as equal to the ensemble averaged concen- 
tration c(x) that we can therefore calculate directly from 
the particle histogram. To test this assumption, c(x) of 
our three-dimensional pores calculated by EFM is shown 
in Fig. 2 by the filled symbols. It can be seen that in- 
deed c(x) agrees very well with the profile calculated by 
relaxing the system into an equilibrium situation. The 
time-consuming relaxation into a stationary state there- 
fore can still be avoided and the fast method of Ref. [12] 
is maintained and improved. 

Examples and results. — In the following, we show 
c(x) and the uncorrected and corrected diffusion coeffi- 
cients D s and D t for several different systems, including 
smooth and rough linear 3d pores (Fig. 4) and different 
Id random walks (Figs. 5 and 6). All diffusion coefficients 
are in units of h 2 /r since we measure the length in units of 
h and the time in units of r. We show that using our new 
method, in all considered cases, both diffusion coefficients 
are in excellent agreement. 

We start with the 3d pores as the experimentally rel- 
evant case (see Fig. 1). We know from [13] that their 
jump lengths I asymptotically show a Levy-distribution 
P{\l/h\) ~ |///i|" (1+/3) with j3 = 3. Therefore we include 
simulations in Id with a similar Levy jump length distribu- 
tion with /3 = 3 into the same figure. To make the simula- 
tions in 3d and Id more comparable we used a composed 
jump length distribution in our Id simulations: Jumps 
larger than h are Levy distributed, jumps smaller than h 
occur equally often. This is to imitate the influence of the 
pore geometry for small jump lengths. Figure 4 shows (a) 
the resulting concentration profiles c(x) for these Id and 
the smooth 3d systems, and the resulting diffusion co- 
efficients (b) without correction and (c) with correction. 
Figures 4(d)-(f) show the equivalent results for rough 3d 
pores of first and third generation. All deviations between 
D s and the uncorrected D t are about 5% and we can see 
they disappear by using EFM. 

It is interesting that the deviations are smaller for the 
rough pores than for v = 0. We believe that this is due to 
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Fig. 4: (a) Particle concentration c(x/L) plotted vs. the nor- 
malized position x/L along the z-axis for smooth 3d pores 
(squares) and the corresponding case of Levy distributed jumps 
in Id (x). (b) Self- (open symbols, +) and transport (filled 
symbols, x) diffusion coefficients for the 3d pores and the Id 
jumps without correction and (c) with correction, (d) Con- 
centration profile and (e) uncorrected and (f) corrected self- 
and transport diffusion coefficients (open and filled symbols, 
respectively) for rough 3d pores of generation v = 1 (triangles) 
and v = 3 (circles). The results were generated from 10 7 runs. 



the very large jump lengths occuring in the smooth pores 
that are stronger suppressed by the boundary roughness. 
Since higher generations only influence smaller length 
scales this effect saturates for higher roughness. Addi- 
tionally for small lengths L and for small times t we see 
huge finite size effects, but both D t and D s are converging 
asymptotical. In this regime the rough geometries show 
more statistical fluctuations since fewer particles diffuse as 
far into the pore. 

Next, we show the calculations of c(x) by EFM for some 
more theoretical cases of Id random walks with differ- 
ent jump length and velocity distributions. In Fig. 5 we 
compare c(x) for these cases as calculated by relaxation 
into a steady state (open symbols) to c(x) as calculated 
by EFM (filled symbols). Also here, the agreement be- 
tween both methods is excellent. In Fig. 5(a) the jump 
lengths I are Gaussian distributed, whereas the velocities 
u arc (i) constant, (ii) proportional to I and (iii) Gaus- 
sian distributed. Figure 5(b) shows c(x) for random walks 
of step size I = nh 7 n £ {1,2,3} with constant velocity 
u = h/r, where u can be identified with the mean velocity 
in real systems. Figure 5(c) shows c(x) for the same ran- 
dom walks but in this case per time step r one complete 
jump is made and jumps of different length therefore have 
different velocities. 

In all considered Id systems, only for systems with a 
single jump length / = ±ft, c(x) docs not deviate from 
the expected linear concentration profile with dc(x) jdx = 
— 1/L. A very pronounced deviation can be seen for sys- 
tems with high probabilities to make either short or long 
jumps per time step. In this case, the concentration pro- 
file can be much higher or much smaller than in the other 
cases (Fig. 5(a) and (c)). For constant velocities or for 
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Fig. 5: Concentration densities c(x/h) of several linear systems 
calculated by relaxation into a stationary state (open symbols) 
and by EFM (filled symbols), plotted versus x for systems of 
length L — 2Ah. (a) Gaussian distributed step sizes with dif- 
ferent velocity profiles: constant velocity (squares), velocity 
proportional to the jump length (circles) and Gaussian dis- 
tributed velocities (triangles), (b) Discrete random walks with 
fixed velocity u = h/r and stepsizes / = lh (squares), \h,2h 
(circles) and lh, 2h, 3h (triangles), (c) Same random walks but 
with velocity u — 1/t. The expected c(x) = Co — x/L only 
applies to the cases of I = lh (in (b) and (c)). The results were 
generated from 10 6 runs. 

a Gaussian velocity profile which is not correlated to the 
jump lengths, the enhancement of the concentration pro- 
file is comparable to the results of the 3d pores. 

As a last example, Fig. 6(a) shows the concentration 
profile for random walks with discrete step sizes of I = lh, 
I E {lh, 2h}, I e {lh, 2h, 3h}, each with equal probability, 
and for Gaussian distibuted step sizes with mean [i = 
and standard deviation a = 2.5h, all with constant u. Fig- 
ure 6(b) shows D s and D t , where D t is calculated using 
Eq. (6) . We see clearly that all cases of non- linear c(x) lead 
to large deviations between both coefficients. Accordingly, 
we now correct D t by using dc(x)/dx from EFM instead 
of Eq. (5). The corrected results are shown in Fig. 6(c). 
For the discrete cases, we also calculated D s and D t ana- 
lytically, yielding D s = D t . These values are displayed in 
Fig. 6(b,c) by straight lines and it can be seen that only 
D t as calculated by the uncorrected procedure deviates 
from the analytical value. 

Finally, we show an easy example to understand the oc- 
curence of a non-linear c{x) close to the boundaries. We 
note that a stationary current requires c(x, t) = const in 
time. Therefore at each lattice site, the number of in- 
and outgoing particles has to be statistically equal and 
constant. For illustration, we have plotted a hypotheti- 
cal linear concentration c(x) = Co — x/L in Fig. 7 for a 
linear system with discrete jumps of lengths I = ±h and 
±2h, each with a probability 1/4. In accordance with the 
method of [12], c(0) is kept constant. For simplicity, we 
assume one jump per time step, as in Fig. 5(c). 

Only for x > 2h, this concentration profile leads to a 
constant current j at all positions x, as a prerequisite for 
stationarity. For the boundary sites, as e.g. for the site at 
x = h, however this is not true: as a consequence of the 
different spectrum of jump probabilities, a strictly linear 
c(x) would not lead to a stationary current at the bound- 
aries. Hence, close to the channel entrance, stationarity 
implies concentration profiles deviating from linearity as 



Fig. 6: (a) Concentration density c(x/L) for random walks with 
discrete jump lengths I = lh (squares), I G {lh, 2h} (circles) 
and I 6 {lh,2h,3h} (triangles up) and Gaussian distributed I 
(triangles down) with /i = and a = 2.5h, plotted versus the 
normalized position x/L. Self- (empty symbols) and transport 
(filled symbols) diffusion coefficients for some of the Id ran- 
dom walks from (a) without correction (b) and with correction 
(c). The straight lines show the values of D s and D t obtained 
analytically. The results were generated from 10 6 runs. 

observed in Fig. 5. Note that the condition of all particles 
starting at x = is different from the given experimental 
situation where particles could enter from outside the sys- 
tem, but it is a necessary ingredient of the method of [12]. 
Even when a homogeneous concentration outside the pore 
is applied a similar behaviour would emerge. 

Conclusions. — We have considered molecular diffu- 
sion in 3d channel pores with different roughness under 
the so-called Knudscn conditions, i.e. for negligible mu- 
tual molecular collisions and for flight times notably ex- 
ceeding the periods of interaction with the pore walls and 
Id systems with many different jump lengths and jump 
time distributions. In the case of the three-dimensional 
pores, which reflect a situation possibly occurring in the 
real nano-world, complete compatibility with the laws of 
normal diffusion is predicted. For non-interacting parti- 
cles, as implied in the considered case of Knudscn diffu- 
sion, this has to lead to equivalence between transport dif- 
fusion and self-diffusion, as illustrated in [2], on the basis 
of Fick's and Einstein's diffusion equations. Assuming a 
constant concentration gradient, errors between ~ 5% (for 
the three-dimensional pores) and ~ 50% (for certain ran- 
dom walks) occur in the calculation of D t when applying 
the method of [12] without calculating c(x) explicitely. 

This complication is related to the fact that, as a typical 
feature of this type of simulations, for molecules entering 
the system only the cross section at x = is considered, 
while jumps out of the system may get to positions far 
beyond this plane. We have shown that the transmission 
probability ft from [12] may as well be applied to calculate 
the real concentration profiles, so that one is released of 
the computational expenses needed for the establishment 
of stationary conditions. With these thus calculated con- 
centration profiles, complete equivalence of both diffusion 
coefficients D s and D t is found in all considered cases. 
Accordingly transport diffusion and self- (or tracer-) dif- 
fusion are equal with respect to both their absolute values 
and their dependence on the surface roughness. Having 
clarified the diffusion behaviour in single pores we have 
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Fig. 7: Sketch of the diffusion process, when jumps of lengths 
±h, ±2h occur with equal probability (one jump per time step) 
and a hypothetical linear profile c(x/h) is applied with fixed 
c(x = 0) = 100. The left part shows the jumps at the left 
boundary, where the postulated linear concentration profile 
would lead to a non-stationary situation. The middle part 
illustrates the flux at a location x = 12h (far from boundary 
influences), where a constant concentration gradient leads to a 
stationary situation. 



and different time intervals, as e.g. changes between two 
successive time steps, but the results remained unchanged. 
[15] Malck K. and Coppens M.-O., J. Chem. Phys., 119 (2003) 
2801. 



established a sound theoretical basis for the exploration of 
mass transfer in the numerous mesopores of nanoporous 
materials. This is important since in many technological 
applications [2,3,5,6], it is this process of mass transfer 
which decides about the performance of these materials. 
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